####################
## Replication file for Bas & Schub "Peaceful Uncertainty"
## Generates Table A11 in Supporting Files
## Cluster Robust Standard Errors (Aronow et al. 2015)

library(foreign)
library(sandwich)
library(Matrix)
rm(list=ls())

## Cluster robust variance estimation function
robust.se.nodfc <- function(model, cluster){
	require(sandwich)
	require(lmtest)
	M <- length(unique(cluster))
	N <- length(cluster)
	K <- model$rank
	dfc <- 1
	uj <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum))
	rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
	rcse.se <- coeftest(model, rcse.cov)
	return(list(rcse.cov, rcse.se))
}


## Set appropriate working directory
triangle2<-read.dta("BasSchub_ISQ_Shocks.dta")

##################
#### Table A11 Model 1
triangle<-triangle2[triangle2$basic==1,]

regSpec <- paste("war2~", 
				paste("Pshock","relcap","contig","s_un_glo",
					"peaceyrs","peaceyrs2","peaceyrs3",
				sep="+"), sep="")

regForm <- as.formula(regSpec)
fit <- glm(regForm, data=triangle, x=T, family="binomial")
summary(fit)

index.sc <- unique(c(triangle$initiator, triangle$target))
index <- index.sc[order(index.sc)]
dyad.mat <- triangle[,c("initiator","target")]
dim(dyad.mat)


# Dyadic cluster robust via multiway decomposition
for(i in 1:length(index)){
	iUp <- index[i]
	clusUp <- apply(dyad.mat, 
					1,
					function(x)as.numeric(iUp %in% x))
	clusIndexUp <- clusUp*-99 + (1-clusUp)*1:nrow(dyad.mat)
	if(i==1){dcrUp <- robust.se.nodfc(fit, clusIndexUp)[[1]]}
	if(i>1){dcrUp <- dcrUp + robust.se.nodfc(fit, clusIndexUp)[[1]]}	
	cat(paste(iUp, " ")); flush.console()
}
# substract naive CR:
dcrUp2 <- dcrUp - robust.se.nodfc(fit, triangle$seclust)[[1]]
# substract HR:
Vhat <- dcrUp2 - (length(index)-2)*vcovHC(fit,type="HC0")

# Cluster Robust SEs and Coefficients
sqrt(diag(Vhat))
summary(fit)$coefficients[,1]


##################
#### Table A11 Model 2
angle<-triangle2[triangle2$interact==1,]

regSpec <- paste("war2~", 
				paste("Pshock","powershift","Pshiftshock","relcap","contig","s_un_glo",
					"peaceyrs","peaceyrs2","peaceyrs3",
				sep="+"), sep="")

regForm <- as.formula(regSpec)
fit <- glm(regForm, data=angle, x=T, family="binomial")
summary(fit)

index.sc <- unique(c(angle$initiator, angle$target))
index <- index.sc[order(index.sc)]
dyad.mat <- angle[,c("initiator","target")]
dim(dyad.mat)


# Dyadic cluster robust via multiway decomposition
for(i in 1:length(index)){
	iUp <- index[i]
	clusUp <- apply(dyad.mat, 
					1,
					function(x)as.numeric(iUp %in% x))
	clusIndexUp <- clusUp*-99 + (1-clusUp)*1:nrow(dyad.mat)
	if(i==1){dcrUp <- robust.se.nodfc(fit, clusIndexUp)[[1]]}
	if(i>1){dcrUp <- dcrUp + robust.se.nodfc(fit, clusIndexUp)[[1]]}	
	cat(paste(iUp, " ")); flush.console()
}
# substract naive CR:
dcrUp2 <- dcrUp - robust.se.nodfc(fit, angle$seclust)[[1]]
# substract HR:
Vhat <- dcrUp2 - (length(index)-2)*vcovHC(fit,type="HC0")

# Cluster Robust SEs and Coefficients
sqrt(diag(Vhat))
summary(fit)$coefficients[,1]

